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We perform N-body simulations of theories with infinite-volume extra dimensions, such as the 
Dvali-Gabadadze-Porrati (DGP) model and its higher-dimensional generalizations, where 4D gravity 
is mediated by massive gravitons. The longitudinal mode of these gravitons mediates an extra scalar 
force, which we model as a density-dependent modification to the Poisson equation. This enhances 
gravitational clustering, particularly on scales that have undergone mild nonlinear processing. While 
the standard non-linear fitting algorithm of Smith et al. overestimates this power enhancement on 
non-linear scales, we present a modified fitting formula that offers a remarkably good fit to our 
power spectra. Due to the uncertainty in galaxy bias, our results are consistent with precision power 
spectrum determinations from galaxy redshift surveys, even for graviton Compton wavelengths as 
small as 300 Mpc. Our model is sufficiently general that we expect it to capture the phenomenology 
of a wide class of related higher-dimensional gravity scenarios. 



I. INTRODUCTION 

One of the most pressing questions in cosmology is 
whether the mounting evidence for dark energy is in fact 
a consequence of a breakdown of Einstein's gravity on 
the largest scales. While observations are converging on 
a background history consistent with that predicted by 
A-Cold Dark Matter (ACDM) cosmology, the most strin- 
gent tests on the standard gravity/ ACDM paradigm will 
come from the formation of large scale structure [1_ . 

In this paper we study structure formation within a 
well-motivated class of infrared- modified gravity theories, 
namely those in which 4D gravity is mediated by massive 
gravitons with characteristic mass r~^. Aside from being 
phenomenologically interesting to study, these theories 
shed new light on the cosmological constant problem [5]. 
Because the graviton is massive, long wavelength sources 
— such as vacuum energy — may effectively decouple 
from gravity, or degravitate [21 SI [S] . 

The simplest and best-studied example in this class 
is the Dvali-Gabadadze-Porrati (DGP) scenario 6 , con- 
sisting of a 3-brane with one extra dimension. Since 
the bulk is fiat and infinite in extent, gravity on the 
brane does not reduce to general relativity at low en- 
ergies. Instead, the force law is approximately l/r^ at 
short distances but weakens to at distances much 
greater than Vc- Much effort has been devoted re- 
cently to confronting DGP with cosmological observa- 
tions [HIHlinilllllllinillllllllSlIIS], most of which 
pertains to the self-accelerated or unstable branch. In 
this work we focus exclusively on the normal or stable 
branch of DGP. 

Our analysis more generally encompasses extensions of 
DGP to higher dimensions — a class of models that are 
non-self accelerating and free of ghost-like instabilities. 
In particular, we are interested in the cascading grav- 
ity framework [TTJ [TSl In this construction, our 3- 
brane lies within a succession of higher-dimensional DGP 
branes, embedded in one another within a flat bulk space- 
time [T71 [TH]. This hierarchy of branes in turn leads to a 
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FIG. 1: Fractional difference between the power spectra for 
degravitation/cascading models (a — 0) with graviton Comp- 
ton wavelength of rc — 300 Mpc (blue dashed) and rc — 500 
Mpc (black solid) and that of standard gravity, without nor- 
malizing to data. The dash-dotted line is the expected differ- 
ence from linear perturbation theory. The dotted line is the 
expected difference assuming the Smith et al. procedure [28| 
for including nonlinear effects. 

force law that successively goes through a 1/r^ regime, 
followed by l/r^, then l/r** etc., as one probes larger dis- 
tances from a source. (A similar cascading behavior was 
also obtained through a different construction in [20;.) 
It was recently argued [3T] that the cosmological pre- 
dictions of cascading gravity models can explain certain 
anomalies in the data. 

When massive, gravitons propagate 5 polarization 
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states: the 2 helicity-2 states of Einstein gravity, 2 
helicity-1 states, and 1 helicity-0 or longitudinal mode. 
The last, traditionally denoted by tt, is most interesting 
phenomenologically — it mediates an extra scalar force 
which enhances structure growth. However, this extra 
force is only effective at sufficiently low density: thanks 
to the Vainshtein effect non-linearities decouple tt 
near astrophysical sources, thereby insuring consistency 



with solar system constraints. (This is qualitatively sim- 
ilar to the chameleon mechanism [23l |24j [25], at play 
in phenomenologically consistent examples of f{R) 
gravity models PT.) Translated to the cosmological con- 
text, this screening mechanism implies that gravity be- 
comes stronger only at late times and on sufficiently large 
scales. 




FIG. 2: Kinetic energy density at z = in a slice of depth 31.25 /i^^Mpc from a pair of 400 ft^^Mpc simulations with the 
same initial conditions, evolved according to standard gravity on the left and a degravitation/cascading model (a = 0, rc — 300 
Mpc) on the right. The units displayed on the scale are arbitrary but common to both panels. We plot kinetic energy density, 
rather than simple overdensity, because the density enhancement in these models is too subtle to detect readily by eye in such 
a plot. Kinetic energy density, however, is enhanced significantly by the greater gravitational force felt by the particles. 



A. Summary of Results 

To perform N-body simulations of DGP and cascad- 
ing/degravitation models, we work in the Newtonian 
limit where the perturbation equations are local on the 
brane. In the DGP context, these have been worked 
out explicitly by Ref. [7] for the case of a spherical top- 
hat perturbation. From the parametric dependence of 
these results, we can infer their generalization to higher- 
dimensional cascading models. For simplicity, we as- 
sume a background expansion history identical to that 
of ACDM cosmology, allowing us to focus on the effects 
from the modified growth history. Moreover, the ACDM 
background cosmology is a realistic approximation for the 
expansion history expected in higher-dimensional DGP 
models [2T], as we review below. 

Figure [l] shows the relative difference in power spec- 



trum between cascading/degravitation models and stan- 
dard gravity, with Vc = 300 (dashed blue line) and 
500 Mpc (black solid line). The spectra from our sim- 
ulations, spanning the range 0.01 ^ k < 1.2 hMpc~^, 
have been extended to smaller k using results from direct 
numerical solution of the linear perturbation equations. 
Because of the extra scalar force mediated by tt, grav- 
itational clustering is enhanced substantially on a wide 
range of scales, and further enhanced on scales that have 
undergone nonlinear processing. And because this ex- 
tra force becomes active at earlier times for smaller rc, 
the enhancement is greater for = 300 than 500 Mpc. 
For comparison, we have also plotted the enhancement 
expected from linear theory (dash-dotted lines) and the 
Smith et al. algorithm [28] to account for non-linear ef- 
fects (dotted hues). 

The enhancement of large scale structure growth in our 
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FIG. 3: Power spectra for rc = 300 Mpc (blue dashed), 
Tc = 500 Mpc (black solid) and standard gravity (red dash- 
dotted), each separately normalized to the Sloan Digital Sky 
Survey main galaxy data _29_ (data points) .This figure demon- 
strates that the apparent difference among the power spectra, 
as viewed via galaxy redshift surveys, is small, given the un- 
certainty in galaxy bias. Compared with the bias of the stan- 
dard gravity power spectrum, the modified spectra biases for 
r-c = 300 (500) Mpc are 57% (65%) relative to that needed 
for the standard gravity results. 



modified gravity simulations is seen directly by making a 
contour plot of the kinetic energy density in degravitated 
and standard gravity simulations. Figure |2] compares the 
results at z = for slices of depth 31.25 ft, Mpc from 400 
Mpc simulations. Structure is clearly more evolved 
in the modified gravity panel due to the 7r-mediated force. 
We plot kinetic energy density, rather than simple over- 
density, because the density enhancement in these mod- 
els is too subtle to detect readily by eye in such a plot. 
Kinetic energy density, however, is greatly enhanced by 
the greater gravitational force felt by the particles in the 
simulation. 

Due to uncertainties in the bias between galaxy red- 
shift surveys and the CDM power spectrum we model, 
however, this enhancement is nevertheless consistent 
with power spectrum determinations such as the Sloan 
Digital Sky Survey (SDSS). Figure [s] shows the fit to 
the Sloan main galaxy power spectrum [29 , assuming 
a scale-independent bias. While the fiducial cosmolog- 
ical parameters assumed here — fim = 0.3, rig = 1.0, 
h = 0.7 — offer a poor fit to the data, even for the stan- 
dard ACDM model, the purpose of this comparison is 
to demonstrate that the apparent difference in P(fc), as 
viewed via galaxy redshift surveys, is small, given the 
uncertainty in galaxy bias. For the record, the bias re- 
quired for our Tc = 300 (500) Mpc models to fit these 
data is 56% (67%) of the bias parameter needed for the 



standard gravity results from our simulations. This is 
a consequence of the greater amplitude of the modified 
gravity CDM power spectra. Note that an important as- 
sumption in this comparison is that of a constant bias. It 
was argued recently that modified gravity theories gener- 
ically lead to a scale-dependent bias [30], an effect we are 
currently quantifying using our simulations. 



II. PHENOMENOLOGY OF 
MASSIVE/RESONANCE GRAVITY 

A. Fundamentals 

The modified gravitational law in standard DGP fol- 
lows from the graviton propagator l/{k'^ + r~^k), where 
Tc sets the cross-over scale from the 4D (l/r^) to the 5D 
(1/r^) regime. While the more complicated nature of 
the force law in higher-dimensional degravitation models 
does not lend itself to such a simple form for the propa- 
gator, a useful parameterization is |5j I31j 



fc2 + r-^(i-")po' ^ ' 

with standard DGP corresponding to a = 1/2. This 
power-law parameterization is not only simple in form, 
but it also makes contact with the far infrared limit of the 
cascading/degravitation propagator: since these models 
all have D > 6 space-time dimensions, the force law 
scales as l/r^~^ in the far infrared, corresponding to a 
propagator that tends to a constant {D > 6) or behaves 
as log k {D = 6) as k ^ 0. Thus all higher-dimensional 
extensions of DGP correspond to a w theories in the 
infrared ^TE . 

In this work, therefore, we shall be primarily inter- 
ested in a = 1/2 and a — 0, corresponding respectively 
to standard DGP and cascading/degravitation models. 
More generally, a is a free parameter within the allowed 
range < a < I. The upper bound is required in order 
for the modification to be relevant in the infrared; the 
lower bound follows from unitarity |31j . 

The above propagator describes a resonance graviton 
— a continuum of massive gravitons — whose spectral 
density peaks at the tiny scale r~^. By virtue of being 
massive spin-2 particles, each graviton state propagates 5 
polarizations, including a helicity-0 or longitudinal mode 
TT. (Cascading gravity models also have extra 4D scalar 
degrees of freedom inherited from the higher-dimensional 
massless graviton [171 132] . For simplicity, in this work we 
shall ignore these scalars and focus exclusively on tt.) 

As emphasized above, this longitudinal mode is re- 
sponsible for nearly all of the interesting phenomenology 
of the models considered here. At the linearized level, 
it contributes an additional T^T'^/6 to the one-particle 
exchange amplitude between conserved sources, which at 
first sight would seem grossly to violate solar system con- 
straints [33] , 
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As Vainshtein [22] realized, however, the weak-field ap- 
proximation is invalid for the longitudinal mode in the 
vicinity of astrophysical sources. Instead, non-linearities 
in TT become important and result in it decoupling from 
matter on scales smaller than a macroscopic scale given 

by m 



l/(l+4(l-a)) 



(2) 



where rsch is the Schwarzschild radius of the source. For 
a = 1/2, this rj,-efFect has been confirmed in explicit 
DGP solutions [Ml US [31] . 

To be more explicit, on scales r <C the leading cor- 
rection to the Newtonian potential is [5J [3T] 



5^) 



r I r 

rsch \rc 



2(l-a) 



(3) 



The above parametric dependence is fixed by two require- 
ments: i) that TT be of order rsch/'^it at r = r^,, to match 
the linear solution; ii) that the solution be analytic in 

Tc for r <C r*. Equation shows that, as desired, 

TT leads to a small correction to the Newtonian potential 
for r <^ Ti,. Thus non-linearities lead to a decoupling of 
TT and recovery of Einstein gravity locally. 

For r ^ r^;, on the other hand, tt yields a correction of 
order unity. 
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(4) 



This is consistent with the tt contribution to the exchange 
amplitude mentioned earlier. 

Thus one can think of the helicity-0 mode as mediat- 
ing a scalar fifth force which is suppressed near high den- 
sity sources but becomes relevant at astrophysically large 
distances. This is precisely the opposite of a fifth force 
mediated by a massive scalar field, where the Yukawa po- 
tential fades away at distances larger than the Compton 
wavelength. The tt behavior instead closely resembles 
that of the f{R) scalar field |27j under the chameleon 
mechanism [2SJ HH [2S] . Indeed we will see that many of 
our results share qualitative features of N-body simula- 
tions of f{R) gravity [37|. 

The key point is that r*, while large for astrophysical 
sources, is nevertheless parametrically smaller than r^. 
Quantitatively, a lO^^M© galaxy has a Vainshtein radius 
of ~ 50 kpc for Tc = 300 Mpc with a = 0. For 
a IO^^Mq galaxy cluster, this gives 1 Mpc. Given these 
scales, we expect tt to play an important role in structure 
formation. 



B. Background Cosmology 

The infrared modifications of gravity discussed above 
should translate in the cosmological context into cor- 



rections to the expansion history. This is certainly the 
case in the standard DGP scenario, where the Friedmann 
equation receives an added contribution proportional to 
H/tc [38'. By virtue of being linear in H, however, this 
correction results in too large a deviation from ACDM ex- 
pansion history, leading to significant tensions with cur- 
rent data fT5]. 

The situation is much more hopeful in higher dimen- 
sional models. As argued in [2T], for general a theories we 
expect power-law corrections to the Friedmann equation 
of the form 



H' = 
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(Recall that in this work we focus exclusively on the 
so-called "normal" branch, as opposed to the self- 
accelerated branch, hence the choice of minus sign on 
the right hand side. The "plus" branch version of this 
equation was introduced in |39j to study generalized self- 
accelerated solutions.) In particular, (|5| agrees with the 
standard DGP Friedmann equation |38j for a — 1/2. 
More importantly, we immediately see that a — (cor- 
responding to 2 or more extra dimensions) leads to an 
expansion history identical to that of ACDM cosmol- 
ogy. Of course, a complete understanding of higher- 
codimension cascading gravity models will inevitably un- 
cover small departures in their expansion history com- 
pared with ACDM, as they do not exactly correspond to 
a = 0. But a — should offer a realistic approximation 
to their modified Friedmann equation since we expect the 
departures to be a slowly- varying function of -ff/r^ |40j . 

In light of these considerations, in this work we there- 
fore assume that the expansion history is identical to that 
of ACDM cosmology. This has the virtue of disentan- 
gling the effects due to a modified growth history, which 
is our primary interest. And while this is only justified 
for higher-dimensional cascading models, we also assume 
a ACDM background expansion for DGP, so that our 
simulations can enlighten us on the a-dependence of the 
modified growth history. 



C. Cosmological Perturbations 

For the purpose of N-body simulations, we are inter- 
ested in sub-Hubble scales and non-relativistic sources. 
In this regime, Lue et al. derived the modified evo- 
lution equations for a spherical top-hat perturbation in 
the standard DGP model; these results were extended 
in |41j . Here we propose a phenomenological extension 
of these results valid for general a theories, which reduce 
to the equations of [7] in the DGP case. 

The change in the perturbation equations can be en- 
tirely encoded as a scale-dependent modification to the 
Poisson equation: 
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(6) 



where e measures the overdensity 6 = Sp/ p in units of rc, 
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and g is defined by 
1 



1 
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where the dot indicates a derivative with respect to 
proper time. Note that the right hand side of Eqn. [6] 
is meant first to be evaluated in real space before being 
transformed into Fourier space to match the form of the 
left hand side. 

A few comments are in order: 

• The above results all agree with [7J for the DGP 
case a = 1/2, with one exception: we have taken 
into account the modified propagator ([T]) in the 
right-hand side of ^ . This mass term has been ne- 
glected in most earlier studies of DGP cosmological 
perturbations [3 121 [13], except in [TS], presumably 
because Vc was assumed to be of order H^^. In this 
work, however, we shall consider much smaller val- 
ues for Tc, and the mass term cannot be dropped. 

• The correction factor on the right hand side of ^ 
encodes the "fifth force" mediated by the longitudi- 
nal mode TT. This is most easily seen by considering 
the linearized approximation, e ^ 1, in which ^ 
simplifies to 
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-AnGpil - g)S . (9) 



The effective correction to Newton's constant, 
Goff — G{1 — g), can be attributed to the behav- 
ior of TT in a cosmological background. Indeed, for 
the universe as a whole both r and rgch can be 
approximated by the Hubble radius, H^^. Thus, 
in the strong coupling regime, Hrc 3> 1, we can 
translate ([s]) to the cosmological context: 
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The parametric dependence agrees precisely 
with (jsj) in the Hrc > 1 limit: g{Hrc > 1) 
l/{H'r^^^^^"\ Similarly, in the weak coupling 
regime, we have g{Hrc ^ 1) = —1/3, which is 
consistent with Q. 



The non-linear dependence on e in ([6]) encodes the 
Vainshtein effect for local overdensities. Indeed, for 
sufficiently large overdensities such that e ^ 1, the 
right-hand side of ^ reduces to its standard grav- 
ity form. The dependence of the effective gravita- 
tional force as a function of 5 is shown in Fig. |4]for 
z = 1. Because these are plotted at z = 1, the en- 
hancement at small S is not yet maximized, thanks 
to the global Vainshtein effect. In other words, g 
has not yet relaxed to its late-time value of —1/3. 
Since g is a slower function of Hvc for a = 1/2 
than for a — 0, the relaxation time is longer for 
a — 1/2, hence the DGP curve has smaller ampli- 
tude at z = 1. 

Figure |4] underscores the fact that our modifica- 
tion depends on the local density, in a way akin 
to the chameleon mechanism in f{R) gravity sim- 
ulations ^7\. This is unlike earlier studies that 
assumed density-independent modifications, such 

as Sana. 
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FIG. 4: EfTective gravitational attraction relative to standard 
gravity at z = 1 as a function of 5p/p. More precisely, what 
is plotted is 1 — 2g{y/l + e — 1 ) /e, which appears on the right- 
hand side of The enhancement compared to standard 
gravity fades away for large Sp/p, a manisfestation of the 
Vainshtein screening effect. 
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• This same non-linear dependence also encodes a 
weak, but interesting, "anti-Vainshtein" effect in 
local underdensities, where the gravitational force 
can become even a bit larger than even its usual 
low-overdensity enhanced value. In other words, if 
Fig. |4] were extended to include underdensities, the 
curves would keep rising slightly in the range S < 0. 
This effect merits further study, but could be rel- 
evant to the void phenomenon, as we will discuss 
briefly in Sec. |V] 

• The analysis of Pfl, on which the above results 
are based, assumes spherical top-hat perturbations. 
For general perturbations, the equation of motion 
for TT in DGP takes a more complicated form and 
is therefore harder to solve. The precise relation 
between our framework and the tt equation of mo- 
tion is spelled out in Appendix]^ Ongoing N-body 
simulations in DGP by R. Scoccimarro have es- 
tablished that a spherical approximation yields a 
power spectrum that agrees well with the full cal- 
culation [44 . 

• Implicit in our modified Poisson equation is some 
coarse-graining for the density field. Otherwise, 
since e is formally a delta function for a point par- 
ticle, our modification would vanish everywhere in 
this extreme case. Of course, since ^ depends 
non-linearly on 6, its solution will unavoidably have 
some sensitivity to the choice of coarse-graining 
scale. As we discuss further in Sec. |IV[ however, 
we will check that this effect is under control by 
comparing simulations of different resolution. Note 
that these averaging issues also pertain to f{R) 
gravity simulations |37| . since the chameleon mech- 
anism depends sensitively on the size and density 
of objects. 

Equations ([6|-([8]) constitute the core ingredients for 
our simulations. Given a density field S and Hubble pa- 
rameter H, we can solve ([6| for the Newtonian potential. 
The latter then specifies how particles evolve in time in 
the usual way. 

We should stress that the above expressions only hold 
well inside the horizon, which is of course all we care 
about for our simulations. In DGP, for instance, the 
corrections to general relativity take on a rather different 
form on super-Hubble scales, as shown by [101 E] and [2] 
in the self-accelerated and normal branch, respectively. 

The generalization of (|9| for near- and super-Hubble 
modes, and the implications for large-scale growth his- 
tory, has been studied recently in |21j . In particular, it 
was shown that modified gravity can explain the lack of 
CMB temperature correlation on > 60 ° angular separa- 
tions |l5l 137] , provided that the curvature perturba- 
tion C is conserved on all relevant scales. The Newtonian 
limit in this case is different than ^ , but the physics is 
qualitatively very similar — enhanced structure growth 
on scales smaller than the graviton Compton wavelength. 



We follow a more conservative approach here. By fo- 
cusing on only two parameters, Tc and a, we are able to 
succinctly explore the large scale structure phenomenol- 
ogy of this class of models. Furthermore, since the 
modifications we introduce are less dependent on the 
somewhat uncertain super-horizon physics of these extra- 
dimensional scenarios, we can expect our general conclu- 
sions to be robust to changes in those scenarios' particu- 
lars. As we learn more from observations and how they 
compare with simulations, we will be able to place better 
constraints on explicit brane constructions. 



D. Linear Regime 

Before launching into the numerical analysis, we can 
gain analytical intuition by studying the linear regime, 
described by (|9|. The physics is illustrated in Fig. 5] 
where the solid lines labeled by aH and rc/a denote the 
Hubble horizon and graviton Compton wavelength, re- 
spectively. 



At early times, H ^ 



gravity is approximately 



standard on all scales since g ~ 0. Once H drops below 
r~^, however, the scalar force mediated by the longitu- 
dinal mode kicks in, and perturbations experience en- 
hanced growth, at least on sufficiently small scales. Thus 
we expect excess power on small scales compared to what 
is expected in ACDM, at least until they reach large over- 
densities. 

On large scales, the graviton mass suppresses the 
growth of modes with k < a/vc- As can be seen from 
Fig. |5] intermediate wavelength modes {k < a/rc) first 
experience a period of enhanced growth from the longi- 
tudinal mode, followed by a period of decay; very long 
wavelength modes (fc <^ a/rc), on the other hand, expe- 
rience only decay from horizon entry until today. 



In a 



Hrc= 1 




Ink 



FIG. 5: Different growth regimes as function of scale factor 
a and comoving wavenumber k. Colored regions correspond 
to modes inside the horizon today. At early times, H > r~^, 
growth proceeds as in general relativity (Green region); for 
H < r^^ and on scales smaller than the graviton Compton 
wavelength {k > a/rc), growth is enhanced thanks to the 
helicity-0 mode (Blue region); on large scales, k > a/rc, the 
graviton mass suppresses growth (Red region). The Vain- 
shtein effect is not included here. 
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III. NUMERICAL SIMULATIONS OF 
MODIFIED GRAVITY 

Our numerical solutions were performed by modifying 
a publicly available Particle Mesh N-body code [48j; a 
modification of the same original code was used in [TB] 
for other modified gravity models. While many more 
accurate methods for computing the formation of large 
scale structure now exist for standard gravity, the particle 
mesh approach is actually better suited to the modifica- 
tion of gravity we consider here. This is because our mod- 
ification, like the particle mesh approximation, does not 
accurately capture the two point interaction between any 
two of the gravitating N-bodies. On the contrary, its im- 
plicit spherical top-hat assumption only works over fairly 
long distance averaging scales, where the coarse-grained 
picture it relies on is a good approximation of reality. 
This is also true of the particle mesh approximation for 
computing N-body evolution. There, too, individual par- 
ticles only source the gravitational potential after being 
coarse-grained over. Recognizing this vindicates our use 
of the less sophisticated particle mesh approach, but will 
prove a difficulty as we attempt to test our model of grav- 
ity with more precision; we will discuss this issue further 
in Sec. |III A[ In short, the nonlinearities in the tt scalar 
field make it both phenomenologically viable and hard to 
simulate. 

The equations derived above are written with an 
implicit assumption of long wavelengths and Fourier 
transforms. This is what allowed us to write the 
d'Alembertian operator as fc^, for instance. Our numer- 
ical solution of the modified Poisson equation (|6| is of 
course performed using discrete Fourier transforms and 
cyclic reduction in the usual way. The effect of the modi- 
fications on the right hand side of ^ is to provide a non- 
linearly transformed density field, with the operations all 
performed in real space. The result of that transforma- 
tion is then converted into Fourier space. The effect of 
the change on the left hand side is even simpler. Since 
we have the Fourier space version of the left hand side, 
we simply use it to solve for <i>fc, albeit using the Green's 
function appropriate for the discrete Fourier transform 
that appears in the code. 



A. The continuum limit and high density screening 

As written, our modified Poisson equation Q does 
not have a smooth continuum limit: as we increase the 
resolution of our simulation, the nonlinearities of the 
equation will cause a larger and larger fraction of the 
mass in the Universe to be screened. Fortunately, how- 
ever, ^ is an approximation to an underlying equation 
of motion which does have a smooth limit. As described 
in Appendix |A] ([6]) is exact for homogeneous density, 
but gives a poor approximation in very sparse regions 
and whenever a significant fraction of the matter lies in 
high-density regions. 



A moment's thought reveals that our modified Poisson 
equation should be applicable provided that the mass in 
typical grid cells are not screened, that is, provided that 
the typical grid smoothing scale, £g, is larger the r^, cor- 
responding to the mass in a given cell. Concretely, let us 
write the density in a given cell as 6p — NM/Iq, where 
M is the mass in a grid cell with mean cosmological den- 
sity, and N is an overdensity factor. It is easily seen that 
our modified Poisson equation is a good approximation 
for cells satisfying 

— - — I I ^ 1- (11) 

go 15 V400 Mpcy ^ ' 

In other words, any grid cell that has an overdensity of 
less than 4x 10^ compared with the average density today 
will not experience substantial screening. This is consis- 
tent with Fig. |4] where the suppression of the extra force 
is significant only for very large overdensities. Note that 
our criterion is independent of the number of particles 
in the simulation because we have averaged over each 
grid cell. This is a good approximation when A^g < 2Np, 
where A'g and A^p denote the number of grid cells and 
particles, respectively, since the Cloud-in-Cell density as- 
signment scheme used in our code smears individual par- 
ticles over several grid cells. 

Of course the above criterion is harder to maintain 
at early times, since the density is larger. But this is 
neither a surprise nor a concern, since extreme overden- 
sities on the length scales of interest are rare until late 
times, cosmologically speaking. Moreover, at early times 
any enhanced force effect is quenched by the cosmological 
screening through the parameter (3. 

Even in the regime where ^ is valid, we can check 
whether a more dense grid leads to a systematic reduction 
in structure growth over the range of scales where we 
trust our simulations. The results of this test, shown in 
Fig.[T4|and discussed in detail in Sec. |IVD] show no such 
systematic effect. 

IV. POWER SPECTRUM 

We perform a variety of runs — approximately 40 repli- 
cates for each parameter combination — each with 512^ 
grid points and 256'^ particles. The model parameter val- 
ues tested in the runs are Tc = 300 and 500 Mpc, with 
a = and 0.5. (The values of Tc quoted here are all in 
physical Mpc, as opposed to h^^Mpc; all other quantities 
throughout the paper are expressed in the usual /i~^Mpc 
units.) We also perform a set of runs with identical initial 
conditions assuming standard gravity. 

Each parameter combination is run in boxes of size 800, 
400, 200, and 100 h^^Mpc. For each of our four simula- 
tion boxes, we can compute a power spectrum. What 
we want, though, is a combined power spectrum that 
includes information over a wide range of scales. The 
power spectrum for a box of a given size L and con- 
taining Np particles, is accurate down to a length scale 
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Power spectrum assembly from 4 boxes 




0.02 0.05 0.10 0.20 0.50 1.00 2.00 
k in h Mpc~^ 



FIG. 6: This plot demonstrates how a full range power spec- 
trum is assembled from our four overlapping simulation boxes. 
The vertical dashed lines indicate the half Nyquist frequencies 
for the four boxes. The orange, black, blue, and red circles 
represent the averaged power spectra for the 800, 400, 200, 
and 100 h~^Mpc boxes, respectively. The simulation pictured 
is that for rc — 500 Mpc in the degravitation model {a — 0). 



2L/Np, or, equivalently, a frequency scale ^Ny = Np/2L 
— the half Nyquist frequency. We therefore join our 
power spectra together in the following way. Starting 
with the largest box, we keep the power spectrum for 
modes between the smallest fc-mode that the box can re- 
solve (e.g., 0.01 hMpc"^ for the 800 h-^Mpc box) and 
the half Nyquist frequency. We then continue the power 
spectrum by taking the next larger frequency mode from 
the next smaller box. That is, the power at fcNy for the 
800 h~^Mpc box is followed in the composite power spec- 
trum by a power value determined from the 400 h~^Mpc 
box. We then keep all of the power spectrum points from 
the 400 h^^Mpc box between that first frequency and the 
400 h^^Mpc half Nyquist frequency. The procedure is re- 
peated to include points from the 200 and 100 /i~^Mpc 
boxes. The power spectrum determined from our sim- 
ulations thus ends at the half Nyquist frequency of the 
100 h~^Mpc box. The constituent parts and the result- 
ing power spectra are illustrated in Fig. |6]for degravita- 
tion model {a — 0) with rc — 500 Mpc. The transition 
between simulation boxes is marked by jumps in scat- 
ter among runs and hence variance. This is a result of 
our decision not to fix the overall normalization of power 
within the smaller boxes, but to let it be randomly sam- 
pled. The modified simulations show larger scatter due 
to their enhanced growth, which exaggerates the differ- 
ences between different initial conditions among small 
box realizations. 



Standard gravity simulations versus Smith et al. 
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FIG. 7: For standard gravity, this shows a comparison of 
A'^{k) = k^P{k)/2TT^ between our power spectrum and that 
obtained through the Smith et al. non-linear fit [28], which 
is represented by the solid black line. The dashed line is the 
linear prediction. 



For standard gravity, we can test the accuracy of our 
simulations by comparing our power spectrum with the 
non- linear fitting algorithm of Smith et al. [28\. The re- 
sult, displayed in Fig. [7j shows very good agreement. 
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Degravitation (a=0.0) DGP (a=0.5) 




0.01 0.02 0.05 0.10 0.20 0.50 1.00 0.01 0.02 0.05 0.10 0.20 0.50 1.00 

k in h Mpc~^ k in h Mpc"^ 



(a) a = (b) a = 0.5 

FIG. 8: The CDM power spectrum from our N-body codes, assembled from the four simulation boxes in the way described in 
the main text. The blue (black) points are for Tc = 300 (500) Mpc. The red points are for standard gravity. The error bars 
represent the variance among numerical runs. Each has been reconstructed by combining results from four simulation boxes, 
and averaged over a number of realizations. The dotted lines show the results from solving the linear perturbation equations; 
the solid (dashed, dash-dotted) lines are the results from the Smith fitting procedure for rc = 500 Mpc (300 Mpc, standard 
gravity). 



Degravitation (a=0.0) DGP (a=0.5) 




(a) a = (b) a = 0.5 

FIG. 9: Plot of kP{k) for our simulations. The convention for curves and data points is identical to that of Fig. [s] 

I 

A. Comparison virith non-linear fitting formulae tity plotted in Fig. [9] is the power spectrum multiplied 

Figures [8] and |9] compare power spectra for DGP and 
degravitation models with standard gravity. The quan- 
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Modified Haiofit vs. Degravitation (a=0.0) 



Modified Haiofit vs. DGP (a=0.5) 
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FIG. 10: Comparison of our modified non-linear fitting algorithm with power spectra from our DGP (left panel) simulations 
and degravitation/cascading (right panel) simulations, each with rc = 300 (blue dashed curve) and 500 Mpc (black solid curve). 
Again, the error bars represent the variance among numerical runs. The modified haiofit parameters are common to both a 
choices and are listed in Table |l] of Appendix [B| 



by the wavenumber, kP(k), to highlight features on non- 
linear scales. As anticipated, the extra scalar attraction 
due to TT enhances gravitational clustering compared to 
ACDM cosmology. Power is increased substantially on a 
wide range of scales, and even more so on scales that 
have undergone nonlinear evolution. The relative en- 
hancement compared to standard gravity is quantified 
in Figs. [T] and [TT] for degravitation and DGP models, re- 
spectively. 

Also plotted in Figs. |8] and [9] are the results of the 
Smith et al. fitting algorithm [28| in each case. We see 
that the fitting approach agrees very well with our simu- 
lations over the entire range of k for standard gravity, and 
for the linear and weakly non-linear scales in the modi- 
fied gravity models. On more strongly nonlinear scales 
{k ^ 0.2 /iMpc~^), however, the Smith et al. algorithm 
heavily overpredicts power relative to what we find in 
our simulations of modified gravity. This is most easily 
seen in the high-fc "bump" in Fig. [9] spanning the range 
0.2 < fc < 1.2 /iMpc"^ The modified gravity "bump" is 
more pronounced than in standard gravity — as can be 
seen in Fig. [3] — but not as large as what the Smith et 
al. algorithm predicts. 

This mismatch on small scales is attributable to the 
Vainshtein mechanism present in our equations, which is 
not captured by the Smith et al. procedure. In modified 
gravity, sufficiently small-scale structures enter the non- 
linear regime early enough that the 7r-mediated force is 
negligible {g k, 0) at that time. Subsequently, as they be- 
come increasingly non-linear, such perturbations simply 



cease to feel any influence of tt. Sufficiently small-scale 
modes thus never experience the effects of modified grav- 
ity. Intermediate-scale perturbations, on the other hand, 
experience a phase of 7r-enhanced growth before going 
non-linear. 

This leads to a scale-dependent pattern of structure 
growth that is different from that expected from general 
relativity. The Smith et al. algorithm, on the other hand, 
was calibrated against simulations of standard gravity 
that were — apart from the horizon scale — genuinely 
scale-free. The algorithm simply scales the non-linear 
power according to the linear prediction and, as a result, 
overpredicts small-scale power for our model. This mis- 
match is also present in simulations of /(i?) gravity |37j . 
in that case due to the chameleon mechanism. 

Remarkably, it is straightforward to recalibrate the 
Smith et al. formulae to capture the distinctive devel- 
opment of nonlinear structure in our simulations. We 
show the results in Fig. 10 To achieve this, we simply 



rescaled most of the Haiofit algorithm's numerical coeffi- 
cients until we found a set that fit our simulations. The 
parameter adjustments necessary are listed in Table [I] 
of Appendix |B] Because the difference between the DGP 
and cascading/degravitation simulations is small over the 
range of our simulations — and expected to grow smaller 
for larger values of — we have provided a single set 
of fitting parameters. For small values of Vf, these two 
models become increasingly different, so the compromise 
parameters fit less well. The better fitting parameters are 
mentioned in the Appendix. A code for calculating these 
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Degravitation (a=0.0) versus DGP (a=0.5) 
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FIG. 11: The same as Fig. [T] for the DGP parameter choice, 
Of = 0.5; Tc = 300 Mpc is represented by the blue dashed line 
and Tc = 500 Mpc by the black solid line. The dash-dotted 
lines are the expected difference from linear perturbation the- 
ory. The dotted lines are the expected difference assuming the 
Smith et al. procedure [2^ for including nonlinear effects. 



FIG. 12: Comparison of the simulations for DGP {a = 0.5) 
and degravitation/cascading models {a — 0.0). The linear 
theory and nonlinearly corrected linear predictions are plotted 
as dash-dotted and dotted lines, respectively, for each value 
of Tc — 300 (500) Mpc in blue (black). The differences are 
small on the scales of the simulations, but become significant 
on larger scales. 



fitting formulae will be made publicly available jiS] . 

This result at first appears to disagree with another 
study, [16 , that found good agreement between Halofit 
and a variety of weak modifications to gravity, includ- 
ing a phenomenological model that mimicked the DGP 
scalar that is also our focus. Their work, however, fo- 
cused on the case where the change in the scalar force 
was explicitly correlated to an alteration of the expan- 
sion history, and where it mediated a repulsive force. 
This constrained their modification of the growth rate 
to be relatively small and to work against, rather than in 
concert with, the Newtonian force. For very large values 
of Vf., our model also recovers standard gravity smoothly, 
and for the large values of Vc considered in [16 , we would 
similarly find much better agreement with the Halofitting 
algorithm. 



B. Dependence upon model parameters 

Another important question is how the growth history 
depends on a. Figure [12] shows the relative difference be- 
tween power spectra for degravitation/cascading (a — 0) 
and DGP {a = 1/2) model. The difference between them 
are only a few percent over the scales probed by the sim- 
ulations and are thus smaller than the variance in the 
simulations, which are suppressed in the plot. This can 



be seen most clearly in the results from the linearized 
perturbation theory results (dash-dotted lines). The one 
clear difference that is evident both in the linearized re- 
sults and in the simulations is that the difference between 
the two values of a is more pronounced for smaller Tc, 
which is sensible. 

The dependence of our results on the graviton Comp- 
ton wavelength is simpler. Smaller Tc leads to more power 
on intermediate scales, but reaches the "cut off' beyond 
which gravity ceases to operate at larger k. As is im- 
plicit in Fig. [3] over the scales we simulate and for the 
relatively small difference in values of Tc we consider, a 
simple linear rescaling is all that is necessary to account 
for the difference in r^. This rescaling can be approxi- 
mated analytically using Eqn. [6] and the equations gov- 
erning linear perturbation growth (see, e.g., [H]). The 
enhancement over standard gravity scales approximately 
as (//oTc)"'^^, so the linear rescaling necessary to go from 
one Tc to another is approximately (rc[l]/rc[2])°-^^. This 
rescaling matches our simulations quite well (see, for in- 
stance. Fig. [l] where a (5/3)°'^^ would put the two lines 
right on top of one another). This rescaling will even- 
tually break down at small k, however, as we can see in 
the linear solution for small k in Fig. 12 It will also 
fail, along with the rest of our approximations, as we go 
to high k where Vainshtein screening becomes systemat- 
ically large. 
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Vainshtein off (lin) versus on (V) for Degravitation (a=0.0) 



Degravitation (a=0.0, black) and DGP (a=0.5, blue) 
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FIG. 13: A clear illustration of the local Vainshtein screen- 
ing effect, due to large overdensities. This plot was made 
using a 200 Mpc/h box for a degravitation/cascading model 
(a = 0), with rc — 300 Mpc plotted as a blue dashed line, 500 
Mpc plotted as a black solid line. The black vertical dashed 
line shows the half Nyquist frequency for this box. The lo- 
cal Vainshtein effect leads to somewhat less development of 
power in the trustworthy portion of the power spectrum, as 
expected. Since the Vainshtein suppression is only effective 
for high overdensites (5 > 100), however, its strongest mani- 
festation is not captured by the simulations performed here. 

Besides further elucidating how our phenomenologi- 
cal model depends on its input parameters, this result 
demonstrates that it would be possible in principle for 
cosmological observations to discriminate between differ- 
ent kinds of higher dimensional gravity theories, since the 
differences between their predictions can be quite pro- 
nounced, especially on very large length scales. 



C. Vainshtein mechanism and large scale structure 

The Vainshtein screening effect in our class of models 
arises in two ways: i) through a global Vainstein effect 
due to the cosmological background density, encoded in 
the function g defined in ([s]); ii) through a local Vain- 
shtein effect from overdensities, encoded in the non-linear 
e-dependence in ([6]) — see Fig. |4] To disentangle these 
two effects, we performed an additional set of simula- 
tions where, instead of using ([6| as the modified Poisson 
equation, we used its linearized form given in ^. The 
latter does not include the Vainshtein screening from lo- 
cal overdensities, only from the cosmological background. 
In other words, this comparison allows us to isolate the 



FIG. 14: This plot demonstrates the dependence of our re- 
sults on box resolution by showing the fractional difference 
between the full resolution runs and a set of simulations done 
with "half resolution in both grid size and particle number 
(ie. a 256^ grid with 128"^ particles). The black circles shows 
the results for a, rc — 300 Mpc cosmology for a cascading / 
degravitation model, the blue circles the results for a DGP 
model. Note that the results for the finer resolution have 
more power on small scales than the coarser grid shows, even 
though we have restricted this plot to those scales when both 
meshes should give valid results. This gives further evidence 
that higher grid resolutions are not causing significant high 
density screening, as argued in §111 A[ 

effects from local Vainshtein screening. 

The results of this comparison in a representative sim- 
ulation box are shown in Fig. [13] The results for each 
simulation box — and the combined power spectra — 
were similar, but focusing on a single box allows us to 
see the changes due to the local Vainshtein effect. 

The main take-away message from this figure is that 
the power is greater without the local Vainshtein effect, 
as expected intuitively. More interestingly, the power 
excess is larger for = 500 Mpc then for 300 Mpc. In- 
deed, for larger non-linearities in tt are more easily 
triggered, hence the local Vainshtein screening is more 
effective. Quantitatively, local screening requires a large 
overdensity in units of rc, that is, e 3> 1. And from ([7|, 
e ^ 1 is easier to achieve for larger rc- (Of course, for 
much larger values of Vc, such that HoVc S> 1, the global 
Vainshtein screening would still be effective today, and 
the local Vainshtein effect would therefore be irrelevant.) 

Going back to Fig. [13] the reason that the difference 
between the two simulations goes to zero past the half 
Nyquist frequency is because the simulations cease to give 
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a valid determination of the power spectrum beyond that 
point. Intuitively, we expect that the reason why these 
highest frequency modes have more power in the linear 
simulation is because of the enhanced kinetic energy of 
the particles in that simulation: since the mesh cannot 
resolve the inter-particle potential on small scales, we 
can see that in simulations with stronger gravity that 
the two particle correlation falls off more quickly than 
in simulations where gravity is weaker on smaller length 
scales. 



D. Sensitivity to averaging 



As discussed in Sees. |TTC] and |ITTX] our modified Pois- 
son equation (|6| requires averaging the density field. 
Since the equation depends non-linearly on 6, one might 
have expected our results to depend on the coarse- 
graining scale. As discussed in §111 A| the relevant scale 
is set by the grid resolution. Therefore, we can quantify 
the dependence on the scale of averaging by comparing 
simulations of different resolution. Figure [T4| shows that 
degrading the resolution of runs by a half, both in grid 
size and in particle number, results in a < 10% change in 
power. Most importantly, the difference is smallest over 
the intermediate length-scales where our simulations are 
most accurate. This demonstrates that the short dis- 
tance, high-density screening discussed in § III A is not 
propagating its effects to these length scales. Indeed, we 
even see that the finer grid — which resolves more high 
density structures in a universe with hierarchical cluster- 
ing — has more power on small scales than the coarse 
grid, whereas the effect of spurious screening would be to 
decrease the power present at higher resolution. 



E. Confronting observations 

The relatively poor fit of our simulation power spectra 
to the SDSS data shown in Fig. [3]should not cause alarm. 
This is due to several factors, most obviously that our 
chosen cosmological parameters were chosen for model 
comparison purposes, not to fit these data. Also, our 
simulations, which contain only dark matter, will neces- 
sarily match real data poorly. 

With regard to data, the most obvious difference be- 
tween the modified and standard models is that the bias 
needed for the modified gravity models to fit these data 
is very different from the typical bias parameter — 56% 
(67%) for 300 (500) Mpc relative to that needed for the 
standard gravity results. This is a natural consequence 
of our model's enhanced gravitational strength; its conse- 
quences bear further study. At the present, the theory of 
bias in modified gravity models has yet to be worked out. 
The study bias in these models is complicated by the fact 
that they introduce scale dependence into the growth fac- 
tor. There are also indications in related work [21j that a 
more complete phenomenology of a massive / resonance 



graviton model may require a more subtle treatment of 
super-horizon evolution. Both of these effects will intro- 
duce scale dependence into the bias. 

A more pressing concern for these models comes from 
recent precision weak lensing measurements |50j . These 
observations observe a power spectrum normalization of 
CTg (f^m/0.25)°-6'^ = 0.837 ± 0.084 on 6l > 85' angular 
scales, in good agreement with the WMAP measure- 
ment [SI]: (Tg = 0.80±0.04. This would appear at first to 
rule out a significant boost in power spectrum amplitude 
from modified gravity. However, lensing data are not so 
easy to interpret in the context of our modified gravity 
model, and do not provide a strong constraint on it as 
yet. 

A simple way to see this is to compare with the results 
in ^43j, which studies the consequences for weak lensing 
of adding a mass term in the Poisson equation. This is 
of course very different from the density-dependent mod- 
ifications considered here. Nevertheless, as can be seen 
from the growth factor shown in their Fig. 1, our class of 
models is mimicked reasonably well by their Yukawa po- 
tential with QfDorc = —0.2 and inverse mass ~ 0.1—1 Mpc 
(similar to our r^,). Using these values, we see that, ac- 
cording to their Fig. 2, this model is as yet in no conflict 
with the weak lensing data. Furthermore, because of 
the bias issues discussed above, the comparison with the 
Sloan Luminous Red Galaxies also shown in their Fig. 2 
should be ignored. 

Another subtlety to keep in mind is that weak lensing 
observations are sensitive to a different gravitational po- 
tential, usually denoted by This quantity is related 
to the Newtonian potential by 
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(12) 



We can understand this result by noting that, unlike non- 
relativistic particles, photons do not couple directly to 
TT. Hence the gravitational potential relevant for their 
motion must be suppressed by a factor of l/l — g relative 
to the Newtonian potential. 

The expression for <f>„ follows from substituting ( 12 1 
into 



$_(fc) = 



1 



1-f (,/m-i) 



V 2(a-l) 



•47rGp-4 . (13) 

The last factor of inGpa'^Sk/k'^ is what one would use in 
standard gravity to translate 6 into a lensing potential. 
Here, however, the two prefactors each act to suppress 
<i>_ relative to the density perturbation: 

• On scales much larger than the graviton Compton 
wavelength, such that fcr^/a <C 1, the first prefactor 
suppresses 
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• On highly non-linear scales, such that e ^ 1, we 
have 
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Since g ^ —1/3 at late times, the local Vainshtein effect 
therefore results in a 4/3 suppression factor for <&_ on 
these scales. 



V. CONCLUSIONS 

Over the next decade observations of the large 
scale structure will subject the ACDM/standard gravity 
paradigm, thus far remarkably successful at accounting 
for cosmological data, to increasingly stringent tests of its 
predictions. In this work we have studied the possibility 
that cosmic acceleration instead stems from a breakdown 
of Einstein gravity at cosmological distances, due to the 
graviton having a small mass. 

Specifically, we have focused on the normal branch of 
the DGP model as well as its extension to higher dimen- 
sions, called cascading gravity. The longitudinal mode of 
the graviton in these models results in stronger gravita- 
tional attraction which, thanks to the Vainshtein screen- 
ing effect, only kicks in at late times and on sufficiently 
large scales. We have focused exclusively on the "nor- 
mal" branch of these theories, which are ghost-free and 
not self-accelerating. (Note that the effects uncovered 
in this work would be exactly the opposite on the self- 
accelerating branch — the longitudinal mode would me- 
diate a repulsive force, thereby impeding the formation 
of structure.) 

We have presented a simple phenomenological pro- 
cedure for calculating the development of large scale 
structure in the standard DGP model and its higher- 
dimensional, cascading gravity cousins. To focus on the 
effects of modifications to the growth history, we have 
assumed a background ACDM expansion history in all 
cases. This should be a good approximation for the 
actual background in cascading models, where the cor- 
rections to the Friedmann equation is expected to be a 
slowly- varying function of Hrc- 

Our N-body simulations confirm the expectation that 
structure is more evolved in modified gravity than in 
ACDM cosmology: 

• Structure grows faster on large scales, leading to 
enhanced clustering for fixed primordial normaliza- 
tion. 

• As a converse, for a power spectrum normalized 
to the power observed today, our model predicts 
systematically less power in the past than would 
be inferred from a standard ACDM evolution. 

• Nonlinear processing reaches longer length scales 
thanks to enhanced gravitational strength, leading 



to an additional enhancement in power on nonlin- 
ear scales relative to standard gravity. 

• Nonetheless, freedom to marginalize over galaxy 
bias allows our modified gravity model to fit the 
power spectrum derived from the Sloan Digital Sky 
Survey main galaxies. 

While our understanding of these simple results are en- 
couraging, they represent only the beginning of our pro- 
gram for testing this class of modifications of gravity. 
Galaxy power spectra, with their bias uncertainty, are 
not the most discriminating measurements for constrain- 
ing or demonstrating the particular effects of these mod- 
els. 

Several observational probes can constrain or confirm 
the phenomenology of our class of models. For instance, 
bulk flows and weak lensing are sensitive directly to col- 
lections of dark matter. Since gravity does not lens 
light with the same potential as it attracts matter in our 
model, though, lensing will not be as useful in constrain- 
ing this model as it may appear at first. Bulk flows, 
though, give a direct measure of the local matter po- 
tential, and are becoming a more robust observational 
tool. Intriguingly, recent results from the compilation of 
peculiar velocity surveys [52^j and systematically uncer- 
tain, but highly suggestive, kinetic Sunyaev-Zel'dovich 
observations [53"| have provided evidence that the pecu- 
liar motions of galaxies on the largest length scales are 
higher than that expected in the standard cosmological 
paradigm. Our model may be able to account for these 
enhanced flows [S3]. 
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APPENDIX A: RELATION TO DECOUPLING 
ARGUMENTS 

In this Appendix we elucidate the relation between our 
modifled Poisson equation and the tt equation of mo- 
tion derived in a decoupling limit of DGP [55l [56j [57] . 
For simplicity, we neglect the effect of the cosmologi- 
cal background and treat the source as embedded in flat 
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Minkowski space. Hence, g —1/3 and p5 p. More- 
over, we focus on scales much smaller than the gravi- 
ton Compton wavelength: kr^ 3> 1. In this approxima- 
tion, (|6| reduces to 
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The contribution from it in the above is the difference 
from the usual Poisson equation, namely 



AiiGp 
9 

8^ 



M-KGrlp 



1 



327rGr2p 

MirGrlp 



27 



1 



27 



1 



(A2) 



In particular, in the limit Gr^p ^ 1, this simplifies to 

AttG 



(A3) 



The factor of 1/3 relative to the standard Poisson equa- 
tion is consistent with the correction in Q. 

Let us compare these results with the tt equation of 
motion of DGP. In a certain decoupling limit of the the- 
ory [55J inSl |57], the longitudinal mode decouples and 
obeys the equation of motion 



3V27r- 



[(V^vr)^ - (V,Vj7r)2] = inGp . (A4) 



In general, because of the V^Vj-tt ter m, we cannot alge- 
braically solve for V^tt, as implied by (A2|. 



1.0 




G p li 



FIG. 15: Comparison of our approximate solution for tt with 
exact profile in the far-field regime. We do not expect this 
mismatch to be very consequential at the level of approxima- 
tion used in this study because most of the dynamics in our 
simulation are driven by regions whose overdensity Gpr^ <^ 1, 
where the approximate and exact solution are most similar. 

However, within a spherical top-hat of radius R, it is 



easily seen that the solution is given by 
3r2 



167-; 



1 



GAnGr^p 
27 



- 1 



for r<R, 



(A5) 



j which agrees precisely with ( A2 1 . After all this is not sur- 

J prising since we obtained ( A2 ) from the spherical top-hat 

(Al) solutions of Lue et al. [7]. At the origin of this agreement 



is the fact that 7r(r) ~ r within the top-hat, and hence 



(A6) 



Therefore, the underlying assumption in ( A2 1 is that (A6) 
holds in general. 

To see what this entails, consider the solution outside 



the top-hat {r > R). From ( A2 1 , we obtain 



27 



approx I 



GM 



3r2 327rGr2 



P 



UirGrlp 



27 



(A7) 

The solution to the exact equation (A4), on the other 
hand, decreases more slowly, since |V7rcxact| ^ l/r^/2 for 
r <ri,, as can be seen from ([3|, and approaches the l/r2 
profile asymptotically: |V7roxact| GMjZr^ for r ^ r*. 
In other words, far away from the source we have 



VTT; 



approx 



^^cxact 



27 



i2'KGrlp 



GA-KGrlp 
27 



- 1 



(A8) 

For small density, Gpr^ <C 1, this ratio approaches unity. 
For large density, on the other hand, our approximation 
tends to underestimate the far- field effect of tt. This ratio 
as a function of Gpr\. 



is plotted in Fig. [151 as a function of Gpr'^. Hence, our 
simulations offer a conservative approximation to the ac- 
tual modifications induced by the scalar-mediated force. 



APPENDIX B: MODIFICATIONS TO 
NON-LINEAR FITTING ALGORITHM 

The changes to the Smith et al. fitting coefficients nec- 
essary to fit our runs are as follows. Any parameters not 
mentioned are left as in the original algorithm. Since we 
have only one set of simulations, we cannot give the full 
dependence of the parameters on the effective index, n, 
or the spectral curvature; these dependencies will, for our 
purposes, be left the same as in the original algorithm. 
The one alteration we make that is not part of the origi- 
nal algorithm is to change the definition of the parameter 
y from y(Smith) = k/k,, to j/(mod) = Qk/ka-', that is, we 
introduce a new parameter, 3. 

The relevant parameters and their modified values are 
listed in Table |l] These parameters have been calibrated 
only for our small sample of runs. A code containing the 
modified Smith algorithm, with updated fitting coeffi- 
cients if we perform a wider variety of runs, will be made 
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available . We provide a single set of parameters be- 
cause both the DGP and Degravitation simulations can 
both be fit reasonably well by them. There is some dis- 
agreement for the smallest values of rc (see Fig. 12 1, 
though. To fit these, it is necessary to use slightly dif- 
ferent parameters for each case. For Degravitation, the 
small Tc fit requires 3 = 1.02, the correction to /3 to be 
2.0, and the correction to /i„ and j/„ to be 0.8. For DGP, 
the small rc fit requires 3 — 1.05, the correction to f3 to 
be 1.9, and the correction to /i„ and to be 0.95. 



logio 


0.84 logio 


logio bn 


logio 6^" + logic 1.1 


logio Cn 


logio c^*<i-h logio 1.05 


logio 


logio Mn" + logic 0-875 


logio '^n 


logio vl'^ + logio 0-875 


Ctn 


0.8 a^M 




1.95/3^"^ 


g 


1.035 



TABLE I: Parameters for the modified non-linear fitting al- 
gorithm. The superscript "Std" indicates the value of the 
parameter in the standard Smith et al. formula. 
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